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Abstract 

This paper presents the non-linear generalization of a previous work on matrix differential models 1 1 1. It focusses on 
the construction of approximate solutions of first-order matrix differential equations Y'(x) = f(x,Y(x)) using matrix- 
cubic splines. An estimation of the approximation error, an algorithm for its implementation and illustrative examples 
for Sylvester and Riccati matrix differential equations are given. 

Keywords and phrases. First order matrix differential equations, Cubic-matrix splines, Sylvester and Riccati differ- 
ential equations. 

1 Introduction 

A great variety of phenomena in physics and engineering can be modelled in the form of matrix-differential equations. 
Although linear matrix-differential equations, whose numerical solutions using cubic matrix splines were presented 
in HI, are valid for a wide range of applications, non-linear equations are also of great interest. This work generalizes 
the approach of [1], providing a novel scheme to numerically solve non-linear differential matrix equations of the 
first-order. Concretely, in this work we will develop a method for the numerical integration of the first order matrix 
differential equation given by 



Y'(x) = f(x,Y(x)) 

Y(a) = Y a 
where Y a ,Y(t) S C rxq , f : [a,b] x C rx « h-> C x «. 



a<x<b, (1.1) 



Different examples of problem (11.11 ) can be found in 13. Numerical schemes to obtain approximate solutions for ( II . lb 
by means of linear multistep methods with constant steps have been devised in [3|. Although there exist a priori error 
bounds for these methods expressed in function of the data problem, these error bounds are given in terms of an expo- 
nential which depends on the integration step h. Therefore, in practice, h will take too small values. Furthermore, these 
methods require some interpolation techniques in order to get a continuous solution [3 1. 



Generalizing the method proposed for the linear case in [ 1 1, here we elaborate an extension using cubic-matrix splines 
in the numerical approximation for the solutions of ( II . lb - In the scalar case, cubic splines were used in [4] for the res- 
olution of ordinary differential equations obtaining approximations that, among other advantages, were of class < ^ )1 in 
the interval [a,b]. These splines are easy to compute and produce an approximation error of only 0(h 4 ). Recently, this 
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1 



method has been used in the resolution of other scalar problems as discussed in [5 1, and even linear matrix problems 
(see HI). The present work extends this powerful scheme to the resolution of matrix problems of the non-linear type 
O- 

This paper is organized as follows. In section [2] we develop the proposed method, whose algorithm is then given in 
Section[3] Finally, in Sections |4]|5] and Impractical examples are presented. 



Throughout this work, we will adopt the notation for norms and matrix cubic splines as in the previous work [ 1 J and 
common in matrix calculus. Following this nomenclature, we define the Kronecker product of A — £ C mx " an( j 
B £ C rxs , denoted by A ® B, as the block matrix 



a \B 



a m \B 



a\ n B 



ClmnB 



The column- vector operator on a matrix A £ C mx " is given by 

A.i 



vec(A) 



A» ;) 



where A.^ 



aik 



IfY — iyij) G C pxq and X = (xij) £ C mx ", then the derivative of a matrix with respect to a matrix is defined by [6 
p.62and81]: 



dY_ 
dX 



dY 
dxu 

dY 

dx m l 



dY 
dx\ n 

dY 



where 



dY 

dx rs 



dy 



Iq 



dyu 

dy p i dy pq 



If X S C m x " , Y e C" x " , Z e C p x q , then the following rule for the derivative of a matrix product with respect to another 
matrix applies [6, p. 84]: 



dXY dX r . . dY 



(1.2) 



where I q and I p denote the identity matrices of dimensions q and p, respectively. If X £ C mx ",Y £ C" xv ,Z £ C pxg , the 
following chain rule [6 p. 88] is valid : 



dZ 
dX 



dx ® h 



dz 



d [vec(y)] 



If A = (a i; ) £ C mx ", the Frobenius norm of A is [ 12 1 given by: 



u\\ F = JZLH z - 

V <=!/= 1 

The following relationship between the 2-norm and Frobenius norm holds [12]: 

l|A|| 2 <||A|| f <^||A|| 2 . 



(1-3) 



(1.4) 



(1-5) 



2 Proposed general method 



Let us consider the problem 
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Y'{x) = f(x,Y(x)) 
Y(a) = Y u 



a<x<b, (2.1) 



where Y a ,Y(t) G C x «, f : [a,/?] x C rxi > ^ C rx «, / S (T), with 

r = {(x,Y) ; a<x<b, YeC rxc i} , (2.2) 

and / fulfills the global Lipschitz's condition 

\\f(x,Yi) - f(x,Y 2 )\\ <L\\Y l -Y 2 \\ , a<x<b,Y u Y 2 eC rx « , (2.3) 

which guarantees the existence and uniqueness of the continuously differentiable solution Y(x) of problem ( 12. U . see 
Qp.99]. 

Let us consider h = (b — d)/n, n being a positive integer, so that the partition of the interval [a,b] is given by 

&\ajb] = { a = x o < x i < ■■■ < x n = b} , X/c = a + kh , k = 0, 1 , . . . , n . (2.4) 

We will construct in each subinterval [a +kh,a + (k+ l)h] a matrix-cubic spline approximating the solution of problem 
( 12. U . For the first interval [a, a + h], we consider that the matrix-cubic spline is defined by 

Si „ (x) = Y(a) + Y'(a)(x-a) + -Y" \a){x- af + i-A (x-a) 3 , (2.5) 

I [«.«+*] 2! 3! 

where the matrix Aq G C rXl? is a parameter to be determined. It is straightforward to check: 



>, (a)=Y(a), S\ (a)=Y'(a)=f(a,Y(a)) 



To fully determine the matrix-cubic spline we still must obtain Y" (a) and Aq. We consider the functions h\ and h 2 
defined by 

hi:[a,b}^>[a,b] h 2 : [a,b] i-> C rxcl 

/ii(jc) = x /12M = Y(x) 

where Y(x) is the theoretical solution of d2.lt . We describe now f(x,Y(x)) as a composition of functions / and (hi,h 2 ), 
that is, let <j> : [a,b] i-» C* 9 be defined by 

0(x) = [fo^.AaJKx) = /(/ixOcJ.Azfjt)) = /(x,r(je)) . 
Thus, ^ is a real variable function of x, and applying theorem 8.9.2 of [6, p. 170] its derivative takes the form: 

D<j> = D(fo(h h h 2 )) = ((D l f)(h u h 2 ))-Dh l + {{D 2 f){h u h 2 ))-Dh 2 , 

where the partial derivatives of /, D\ (/), D 2 (f) exist and are continuous since it is assumed that / S c € l (T). By ( 12.1b 
it is clear that 

rf < V "7< J » r = [ V ecf { x,Y(x))f . 
ax 

Next, applying the chain rule for matrix functions ( 11.2b and then taking the derivative of a matrix with respect to a 
matrix, ( 11.31 , one obtains 



d/(*,y(*)) , r r vl U1 r^ f i d/(*,y(*)) 



T"W = ^ + [vecf(x,Y(xj)Y ®7, 



<9 vec y (x) 



(2.6) 



We are now in the position to evaluate y"(a) using 
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By imposing that ( 12.5b is a solution of problem ( 12.11 ) in x = a + h, we have: 

S\ (a + h) = f(a + h,S\ (a + h)) , 

\[a,a+h] \ \\a,a+h\ J 



and obtain from (12.71 the matrix equation with only one unknown matrix Aq 
A 



2 



/( a + h,Y(a)+Y'(a)h + ^Y"(a)h 2 + ^A h 3 ) - Y' (a) - Y" (a)h 



(2.7) 



(2.8) 



Assuming that the matrix equation ( 12.81 ) has only one solution Aq, the matrix-cubic spline is totally determined in the 
interval [a,a + h]. 

Now, in the interval [a + h, a + 2h], the matrix-cubic spline takes the form 



(x) = Si (a+h)+S\ (a+h)(x-(a + h)) 

+ -S'{ (a+h)(x-(a+h)) 2 + ^-Ai(x-(a+h)) 3 , 

1 ! I [a,o+A] 3 ! 



(2.9) 



sothatS(jc) is ofclass ^{[a^b]) on [0,0 + /;] U [a +h, a +2h], and all coefficients of the matrix-cubic spline 5| ^ ^ (x) 

are determined with the exception of A\ G C' x ' 7 . By construction, matrix-cubic spline (12.91 satisfies the differential 
equation ( 12. U in x = a + h. We can obtain A \ by requiring that the differential equation ( 12. U holds at point x = a + 2lr. 

S\ (a + 2h) = f(a + 2h,S\ (a + 2h) 

\{a+h,a+2h] V \[a+k,a+2k] 

Expanding, we obtain the matrix equation with only one unknown matrix A\: 



2 
JP- 



f[a + 2h,S\ (a + h)+S\ (a + h)h + -S\' (a + h)h 2 +-A l h 3, 

\[a,a+h} \[a,a+h] 2 \[a,a+h] 6 



S\ (a + h)-S'( (a + h)h 

\[a,a+h] \[a.a+h] 



(2.10) 



Let us assume that the matrix equation J2. 10b has only one solution Aj. This way the spline is totally determined in the 
interval [a + h,a + 2h]. 

Iterating this process, let us construct the matrix-cubic spline taking [a + (k— l)h,a + kh] as the last subinterval. For 
the next subinterval [a + kh, a + (k + 1 )h], we define the corresponding matrix-cubic spline as 



(2.11) 



where 



PkW = Ett^ (a+kh)(x-(a + kh)) k . 

"die! \[a+(k-\)h,a+kh] 



(2.12) 



With this definition, the matrix-cubic spline is S(x) S c € 2 M [a + jh, a + (j + 1 )h] and fulfills the differential equa- 

V=0 / 

tion d2.U at point x = a + kh. As an additional requirement, we assume that S(x) satisfies the differential equation ( 12. j} 
at the point x = a + (fe + l)/z: 

S\ (a+(k+l)h) = f(a+(k+l)h,S\ (a+(k+l)h) 

\[a+kh,a+{k+l)h} \ |[o+tt,o+(*+l)»] v v ' y 
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and expanding this equation with the unknown matrix A k yields 

2 

W 
1 



At 



f I a+(k+ \)h,S\ (a + kh)+S\ (a + kh)h 

* \[a+{k-l)h,a+kh] \[a+(k-l)h,a+kh] 



-S\' (a + kh)h 2 + -A k hr 

2 | [«+(*-!)*,«+**] 6 



S( 



|[n+(*-l)/!,a+*/i] 



(a + kh) — Si (a + kh)h 

|[«+(t-l)A,o+tt] 



(2.13) 



Note that this matrix equation (12.13b is analogous to equations ( 12.81 l and ( 12.10b . when k = and k = 1, respectively. 
We will show that these equations have an unique solution using a fixed-point argument. 



For a fixed /i, we will consider the matrix function of matrix variable g : C x<? i— > C rx<? defined by 



2 

1 



f[a+(k+\)h,S\ (a + kh)+S\ (a + kh)h 

' 1 |[«+(*-l)*,o+*ft] |[ a +(i-l)A,a+f»] 



2 |[o+(*-l)A,o+»] 



(a + kh)h 2 



1 



-ttt 



(2.14) 



- S\ (a + kh)-S\' (a + kh)h 

I [a+(t- 1 )h,a+kh] | [«+(*- 1 )h,a+kh] 

Relation ( 12.131 holds if and only if A k = g^A^), that is, if A k is a fixed point for function g(T). 

Observe that by using ( 12.121 ) and applying the global Lipschitz's condition ( 12.31 ) it follows that 

Lh 

||s(7i)-s(72)||< — ||7i-T 2 ||. 

Taking h < 3/L, g(T) yields a contractive matrix function, which guarantees that equation (12.131 ) has unique solutions 
A k for k = 0,1,..., 7i— 1. Hence, the matrix-cubic spline is completely determined. Taking into account JU Theorem 
5], the following result can be established. 



Theorem 2.1 Let be L the Lipschitz constant defined by ( 12. 3D . If h < 3/L, then the matrix-cubic spline S[x) exists 
in each subinterval [a + kh, a + (k + l)h\ k = 0, 1 , . . . ,n — 1, as defined in the previous construction. Furthermore, if 
f £ ^ 3 (r), then \\Y(x) — S(x)\\ — 0(/j 4 ) Vx £ [a,b], where Y(x) is the theoretical solution of '( 12. 7 P . 



3 Algorithm 

The following algorithm is designed to compute the approximate solution of ( 12.11 ) by means of matrix-cubic splines in 
the interval [a,b] with an error of the order 0(h 4 ) under conditions of theorem |2~T1 

• Determine the constant Y"(a) given by ( 12.61 ). Take n > L(b — a)/3,h=(b — a) jn and 
the partition A[ a & i defined by Eq. ( 12.41 ). 

• Solve the matrix equation ( 12.81 ) for k = and determine St (x) of Eq. ( 12. 51 ). 

• Solve the matrix equation (|2.13t iteratively for k=l,...,n — 1, and then compute the 
splines Si (x) according to Eq. (12. lit . 

| [a+kh,a+(k+l)h] 

Depending on the function f(t,Y), matrix equations ( 12.8b and ( 12.131 ) can be solved explicitly (see |8]) or using the 
iterative method (see for example 13): 

Tf + j = giTf) , where Tq is an arbitrary matrix in C' XI! , s = 0, 1, . . . ,n — 1 
and g(r) is given for ( 12.14b . In the following section, we will test the algorithm proposed. 
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4 Example: A non-linear vector system 

We consider the next non-linear vector differential system 

y\{x) = - 1 + e x - sin (x) + sin (y 2 (x) ) 



y' 2 {> 



i _ i 

4+yi (x) 2 S+e^+le* cos (x)-sin 2 (x) 



> 0<jc< 1 



(4.1) 



yx(0) = 2, y 2 (0) = f 

It is easy to check that this problem has the exact solution y i (x) = e x + cos (x) , y 2 (x) — f , so in this particular case we 
will be able to obtain the exact error of our numerical estimates. 

We can rewrite ( 14.1b in the compact form 

Y'(x) = F(x,Y) 

2 \ \ 0<x<l,Y(x)=( )eF:.f!.v.n 



Y(0) 



thusF'(0)=F(0,( s 



— 1 + e* — sin (jc) + sin (y 2 (x) ) 
l l 

4+J'i (x) 2 5+e 2x +2e v cos (x)-sin 2 (x) 



^ ] . We calculate F"(0) using ( 12.6b . in this case, one gets 



vec(y W )=F W =(^j), 



yi (x) \ dF(x,Y(x)) 



dx 



/ e x — cos (jt) 



2^ +2^ cos (x)— 2e*sin (x) — 2cos(x) sin(x) 



On the other hand, we have 



\ (5+e 2t +2e t cos(.v)-sin 2 (x)) T " 

[vecF(x,Y(x))] T ®I 2 



( — 1 + — sin (x) + sin (y 2 (x) ) -. — Vtt 1 ^~r\ ttt i £3 h 

\ w WA yy 4+ yi (x) 2 5+e 2r +2e r cos(x)-sin 2 (x) ,1 ^ 12 

( (-l+e T -sin(x) + sin(y 2 (x)))/ 2 | (j^- 



5+e 2j[ +2e A 'cos (x)-sin 2 (x) ' ^ 2 



(4.2) 



(4.3) 



— \+e x — sin (x)+sin (>'2(x)) 

-l+r r -sin(x)+sin(y 2 (x)) 



4+yi (J") 2 5+e-->+2c' cos(j)-sin 2 (x) 

! 



and 



5f (x,F(x)) 
<9 vec T(x) 



, dF(x.Y(x)) , 
\ dy 2 I 



( ^-(-l+^-sin(x) + sin(y 2 (x))) \ 



dy\ \ 4+_vi (x) 2 5+g 2l +2e'' cos (xj-sin 2 (x) 

- ( - 1 + e x - sin (x) + sin (y 2 (x))) 
_( 1 1 

\ V4+yi(x) 2 5+e 2j '+2e-' [ cos (x)-sin 2 (x) ) / 



4+vi {x) 2 5+e 2 *+2e' cos (*) -sin 2 (x) 



( \ 

-2>'i(x) 

(4+ViW 2 ) 2 

cos(y 2 (x)) 

V J 



Therefore 

[vecF(x,Y(x))] T ' ®I 2 
and by d4.3b - d4.61 l one concludes 



dF(x,Y(x)) _ ( (4+^x) 2 " 9+2^+4^ cos (x)+cos(2x) ) COS O^M) 



<9 vec 



2 yi (x) ( - 1 +g' - sin (x) +sin (y 2 (x) ) ) 
(4+>'i(x) 2 



(4.4) 



(4.5) 



(4.6) 
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Tntf j r\7cil 
111 LCI Veil 


Approximation 


IVltlA. HjIIUI 


[0,0.1] 


/ 2 + x + 0.177917x 3 \ 
^ f - 5.62424 x 10~ 6 x 3 J 


2.83337 x 10~ 6 


[0.1,0.2] 


( 


i qqqck j_ i nm^Si- n m -iqiaiy^ _i_ n ooAni i 

1 .yyyyj ~r 1 .UU1 jOI — U.U1 jo J4Z.Y -)- \J.ZZH\Jj LX 

1.5708 + 6.67857 x 10~ 7 x- 6.67857 x lO^x 2 +0.0000 166377x 3 


) 


2.83337 x 10~ 6 


[0.2,0.3] 




/ 1 .yyy I J -f- 1 — U.UZV loZZ_Y -f- U.Z4V01 IX \ 

^ 1 .5708 - 4.57386 x 10~ 6 x + 0.00001953x 2 - 0.0000270433x 3 J 




2.94712 x 10~ 6 


[0.3,0.4] 


( 


1 .yyo4i -+- 1 .Ul IojX — U.U / j IvyjlX -f- u.zyyi4Zx 
v 1. 57079 + 0.0000126873x-0.0000380073x 2 +0.000036887 lx 3 , 


) 


2.94712 x 10~ 6 


[0.4,0.5] 




} 1.99655+ 1.0318x-0.108685x 2 +0.328246x 3 S 
v 1.5708-0.0000271633x + 0.0000616192x 2 -0.0000461351x 3 j 




3.0698 x 10~ 6 


[0.5,0.6] 


( 


1 .9899 + 1 .07 17 lx - 0. 188509x 2 + 0.38 1462x 3 
v 1. 57079 + 0.0000485499x-0.0000898071x 2 + 0.0000548 159x 3 , 


) 


3.0698 x 10~ 6 


[0.6,0.7] 




( 1.98277 + 1.10736x-0.247933x 2 +0.414475x 3 \ 
v 1.57081 - 0.0000785694x + 0.000122058x 2 -0.0000628872x 3 j 




3.20977 x 10~ 6 


[0.7,0.8] 




f 1.96307+ 1.19177x-0.368517x 2 +0.471896x 3 \ 
\ 1.57077 + 0.0001 17333x - 0.000157802x 2 +0.0000703796x 3 J 




3.20977 x 10~ 6 


[0.8,0.9] 


/ 1.94382+ 1.26395x- 0.45874x 2 +0.509489x 3 \ 
V 1.57084 - 0.000166 lx + 0.00019649x 2 -0.0000772419x 3 ) 


3.37764 x 10~ 6 


[0.9,1.0] 




f 1.89829+ 1.41574x-0.627395x 2 +0.571954x 3 \ 
^ 1. 57073 + 0.000224533x-0.000237548x 2 +0.0000835 127x 3 J 




3.37764 x 10~ 6 



Table 1: Approximation for vector differential system (14. j} in the interval [0, 1] with step size h = 0.1. 



Y"(x) = 



dF(x 7 Y(xj) 
dx 

e x — cos (x) 



[vecF(x,Y(x))Y ®I 2 
l 



dF(x,Y(x)) 
d vec Y(x) 



4+yi(x) 2 9+2e 2l +4e J: cos(A:)+cos(2j:) 



cos(y 2 (x)) 



2e- x +2e x cos (x) -2e x sin (a)-2 cos (x) sin (x) 2y\ (x) ( - 1 +g r - sin (x) +sin (y 2 (x) ) ) 



(5+c Zr +2e l "cos(x)-sin 2 (x)) r (4+yi (x) 2 )* 

Taking into account that yi (0) = 2,^2(0) = j and evaluating Y"(x) of ( 14. 7t when x — 0, one gets Y"(0) 
It is straightforward to show that F, defined by ( 14.21 ), fulfills the global Lipschitz's condition 

\\f(x,Y) - f(x,Z)\\ < \\Y-Z\\ , 0<x<l ,F,ZeK 2 , 



(4.7) 



(4.8) 



thus, we can take L given by ( 12.31 l as L = 1. Therefore, we need to take h < 3/L and thus h = 0.1 for example. The 
results are generated with Mathematica using FindRoot function to solve the emerging algebraic equations, and are 
summarized in Table 1 . In each interval, we evaluated the difference between the estimates of our numerical approach 
and the exact solution, and then take the Frobenius norm of this difference. The maximum of these errors are indicated 
in the third column for each subinterval. 



5 Example: Sylvester matrix differential equation 

Linear matrix differential equations of the type 
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3.5e-06 



3e-06 




approx. interval 

Figure 1: Representing the Frobenius error margins for vector differential system ( 14. U in the interval [0, 1] with step 
size h = 0.1. 



Y'(x) = A(x)Y(x)+Y(x)B(x)+C{x) 
Y(a) = Y a 



a<x<b, Y(x),A(x),B(x),C(x) E C rxr , 



(5.1) 



arise in many fields of science and engineering. In the case of constant coefficients has been studied by several authors 
(see for example iflOl ). However, the variable-coefficient case has so far received little numerical treatment in the 
literature. We can observe that the proposed method require the matrix functions A(x),B(x) and C(x) to be differen- 
tiable, while, for example, in the method proposed in ifTTl . it is necessary that A(x),B(x) have continuous second-order 
derivatives and C(x) continuous in the domain a <x <b. 



As an example, here let us consider the Sylvester problem ( 15.11 ) with 

= (°. *rv B(x)=(° n *v c( X ) 



F(0) 



1 
1 



Y(x) e C 2x2 , 0<x< 1 



-e- x (l+x 2 ) -2e- x x 
1 — e~ x x —x 2 



(5.2) 



This problem has an exact solution Y(x) 
error of our numerical estimates. 



, so in this particular case we will be able to obtain the exact 



As we have max 

xe[0,i 



xe 
x 



x 




< 1 .69443, one can take the constant L given for ( |2.3t as L = 2. 
Taking the derivative of Y'(x) —A(x)Y(x) + Y(x)B(x) + C(x), gives: 
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Interval 


Approximation 


Max. Error 


[0,0.1] 


/ l-j + 0.5x 2 -0.1612x- j \ 
V x 1 / 


1.33472 x 10~ 6 


[0.1,0.2] 


/ l-0.9994x + 0.4938x 2 -0.1406x J \ 
V x 1 / 


1.33472 x 10~ 6 


[0.2,0.3] 


/ l-0.9984x + 0.4890x 2 -0.1325x J \ 
V x 1 / 


1.2445 x 10~ 6 


[0.3,0.4] 


/ 0.9994 -0.9936x + 0.4728x 2 -0.1 146x J \ 
V x 1 / 


1.2445 x 10~ 6 


[0.4,0.5] 


/ 0.9991 -0.9909x + 0.466 lx 2 -0.1090x J -O.OOOlx 2 \ 
\ i non 1 v o non 1 r z i 1 

\ 1.UUU1X U.UUU1X 1 1 


1.17402 x 10~ 6 


[0.5,0.6] 


( 0.9971 -0.9791x + 0.4426x 2 -0.0933x J -O.OOOlx + 0.0002x 2 -0.0001x J \ 
I n oooor -i- n nnn?v 2 nnoniv-' noooo / 

\ U.7777A t U.UUUZA U.UvUIa VJ.yyyy 1 


1.17402 x 10~ 6 


[0.6,0.7] 


( 0.9963 - 0.9732x + 0.4361x 2 -0.0898x J 0.0002x - 0.0004X 2 + 0.0002x J \ 
\ 1 000? T — O 0004 y 2 +0 000? r 3 1 / 

V A . UUUiA \J . WWA \j . UV7UZ.A 1 1 


1.12331 x 10~ 6 


[0.7,0.8] 


( 0.9916 - 0.9549x + 0.4071X 2 -0.0759br 1 0.0001 -0.0004x + 0.0006x 2 -0.0003x J \ 
^ 0.0001 +0.9996x + 0.0006x 2 -0.0003* 3 0.9999 J 


1.12331 x 10~ 6 


[0.8,0.9] 


( 0.9906 - 0.95 12x + 0.4025x 2 -0.0739x ;i 0.0002 + 0.0007x - 0.0009x 2 +0.0003x J \ 
I -0.0002+ 1.0007x-0.0009x 2 +0.0003x 3 l+0.0001x 2 ) 


1.09412 x 10~ 6 


[0.9,1.0] 


( 0.9816-0.9212x + 0.3691x 2 -0.0616X 1 0.0004 - 0.001 lx + 0.0012x 2 -0.0004x J \ 
^ 0.0004 + 0.0.9989x + 0.0012x 2 -0.0004x 3 0.9999 + 0.0002x - 0.0002x 2 ) 


1.09412 x 10~ 6 



Table 2: Approximation for Sylvester matrix differential equation d5.2t in the interval [0, 1] with step size /z = 0.1. 



Y"(x) = (A'(x) + (A(x)) 2 )y(x)+y(x)((B(x)) 2 +B'(x)^ 

+ 2A(x)Y(x)B(x)+A(x)C(x)+C(x)B(x)+C'(x) . (5.3) 

WeseethatT'(O) = (~i q j , and by applying (O it is F"(0) = ( * jj 

In this numerical example, we take n = 10 such that n > L{b — a) /3 and h = 0. 1 = (b — a)/n. The results are generated 
with Mathematica using the Bartels-Stewart algorithm (see for example [ 12 1) to solve the emerging algebraic equations, 
and are summarized in Table 2, where the numerical estimates have been rounded to the fourth relevant digit. In each 
interval, we evaluated the difference between the estimates of our numerical approach and the exact solution, and then 
take the Frobenius norm of this difference. The maximum of these errors are indicated in the third column for each 
subinterval. 



6 Example: Riccati matrix differential equation 

Rectangular non-symmetric Riccati matrix-differential equation of the type 

Y'(x) = C(x)-D(x)Y(x)-Y(x)A(x)-Y(x)B(x)Y(x) ] 

> 0<x<c, (6.1) 
Y(0) = Y J 

where the unknown Y(x) G C pxc > and coefficients A(f) e C« x «,B(x) E C xp ,C(x) e *«,£>(*) € C pxp are differen- 
tiable matrix-valued functions arise frequently in important applications to classical control theory [13] and as decou- 
pling techniques for both the analytic and numerical study of boundary value problems lfl4l . The Riccati equation ( 16.1b 
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1.4e-06 - 



1.2e-06 




Figure 2: Representing the absolute error margins for the Sylvester matrix differential equation ( 15.21 ) in the interval 
[0, 1] with step size /i = 0.1. 



has been studied extensively, and different resolution techniques have been introduced (see [ 15] and references therein). 



The study of the Riccati equation (16. U is closely related to the underlying linear system 

X'(x) = S{x)X(x) 



X(0) 



l 1 
Y 



where X(x) = 



U(x) 
V(x) 



S(x) = 



A(x) B(x) 
C(x) -D(x) 



Specifying the solution of (16. U is given by 

Y(x) = V(x)U-\x) 
where Y(x) is defined in the interval where U(x) is invertible, see lfl6l . 



(6.2) 



(6.3) 



Taking into account lemma 1 and 2 of (TTJ, U(x) is invertible in the interval [0,5] and the solution Y(x) of problem 
( 16. II ) satisfies 



||F(x)|| < M ,M = (l-5^ exp(5^o)wo) w exp(<5£ ), 
where 8 is a positive number satisfying 

Sk Q + log (8) < -log(<7owo) , 

and 



(6.4) 



(6.5) 



ko = max 



A(x) B{x) 
C{x) -D(x) 



0<x<c 



Wq 



max{||A(;c) B(x)\\ ;0<x<c} 



Y 



(6.6) 
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In accordance with [18. p. 1064], we consider the matrix-valued function 



F(x,Y) = C(x) -D(x)Y -YA(x) -YB(x)Y 



(6.7) 



then, if we define 



and IIYII < M, 



a = mp{\\A(x)\\ ;0 < x < 8} 

b = sup{||B(x)||;0<x<5} 

c = sup{||C(x)|| ;0 < x < 8} 

d = sup{||Z)(jt)||;0<x< 5} 

< M, with M gives by ( 16.4b . the following local Lipschitz condition holds 

= a + d + 2bM . 



(6.8) 





Y-Y 







In addition, if | |F || <N, 



\\F(x,Y)\\ <c+N(a + d + bN) . 



(6.9) 



(6.10) 



Using the proposed spline method, the only one solution of the matrix equations (12.8b and ( 12.13b for k = 1 — 1 
is guaranteed using a fixed-point argument and the global Lipschitz's condition (12.3b . In our case, we need to prove 
the only one solution of the matrix equations ( 12.8b and (12.13b using a fixed point argument and the local Lipschitz's 
condition ( 16.9b . 

We start with the matrix equation ( 12.8b . Let us suppose that ||T*|| <N\. Taking into account (16.10b . we take 

' N 2 = \\Y(a)\\+h\\Y'(a)\\ + ^\\Y"(a)\\ + ^N l 



< JV 3 = 
. N 4 = 



c + N 2 (a + d + bN 2 ) 

|y'( a )||+A||y"(fl)||) 



(6.11) 



witha,b,c given by ([6J]>, and let be N = max { N\ , N 2 , N 3 , N 4 , M} with M gives by ([63]) . Let be si = {Y G C rx q ; \ \ Y \ \ < N } 
and we consider the continuous matrix-valued function of matrix variable g : C rx<? C rx <? defined by ( 12.141 ) for k = 0. 
It is simple to verify that if T G st ', by ( 16.1 lb and ( 16.10b then g(T) £ si . Thus, g : si i— > si and Aq is a fixed point of 
g. In addition, if Ti , T 2 G si , ||7i || < M, \\T 2 1| < M, has then that for / defined by (EZK / fulfills the local Lipschitz's 
condition ( 16.91 ) and taking h < 3/L, g(T) yields a contractive matrix function, which guarantees that equation ( 12.81 ) has 
unique solutions Aq. Hence, the matrix-cubic spline is completely determined in [a,a + h]. 



For k— 1, ... ,n — 1, fixed, supposed construct cubic-matrix spline S(x) taking [a+ (k— l)h,a + kh] as the last subin- 
terval, for the next subinterval [a + kh,a + (k+ 1 )h], to define the corresponding spline we need determine A* G C rx<? 
as the only one solution of the matrix equation ( 12.131 ). Let us suppose that ||T|| <N\. Taking into account ( |6.10b , we 
take 



N 2 = 
N3 = 
N~ 4 = 



S\ (a+kh) 

\[a+(k-l)h,a+kh\ v ' 



S\ (a- 

| [«+(*-!)*,«+**] 



kh) 



S'f (a+kh) 

|[B+(Jk-l)*,0+**] 



4^ 



c + N 2 (a- 
N3 + 



d + bN 2 ) 



(6.12) 



2 

77 



S\ (a+kh) 

\[a+(k-\)h,a+kh\ 



S'f (a+kh) 

|[a+(Jfc-l)*,«+**] 



with a,b,c given by (ISTSl and let be N = max , N 2 , N 3 , N 4 ,m} with M gives by (El). Let be si = jy G C rxq ; \ \ Y 1 1 < j 

and we consider the continuous matrix-valued function of matrix variable g : C rxq i— > C XCI defined by ( 12.141) . 

It is simple to verify that if T G si, by ( 16.12b and (16.10b then g(T) G si. Thus, g : si t—> si and A^ is a fixed point of 

g. In addition, if T x , T 2 € si , ||Ti|| <M,||T 2 || <M, has then that for / defined by do\7l), f fulfills the local Lipschitz's 
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1 n f^i"\/Q 1 
111 LCI Vd.1 


Approximation 


IVlclA. J_/1IUI 


[0,0.01] 


/ l+x + 0.5x 2 +0. 167224 x J \ 

U 2 * ) 


1.39903 x 10~ 10 


[0.01,0.02] 


( 


' l.+x + 0.499933 x 2 +0.169461x J N 

v x2 x 


) 


1.39903 x 10- 10 


[0.02,0.03] 




/ l+x + 0.499864x 2 +0. 17061 x J \ 




1.41977 x 10~ 10 


[0.03,0.04] 




f 1 + 1. 00001 x + 0.49966x 2 +0.172877x 

V x 2 x 


') 




1.41977 x 10~ 10 


[0.04,0.05] 


( 


' 


1 + 1.00001x + 0.499518x 2 +0.174063x ;S N 
x , 


1 


1.44084 x 10- 10 


[0.05,0.06] 


( 


' 


l + 1.00003x + 0.499173x 2 +0.176362X" 1 

X j 


I 


1.44084 x 10" 10 


[0.06,0.07] 


( 


' 0.999999+1. 00004 x + 0.498952x 2 +0.177587x J N 

v x 2 X J 


) 


1.46223 x 10~ 10 


[0.07,0.08] 


( 


' 0.999998 + 1.00008x + 0.498463x 2 +0.179918x J N 

\ x x ; 


) 


1.46223 x 10~ 10 


[0.08,0.09] 






0.999998+1. 0001 x + 0.498 16x 2 +0.181181x j \ 




1.48391 x 10- 10 


[0.09,0.1] 


( 


' 0.999996+ 1.00016x + 0.497521x 2 +0.183546x J N 

, X 2 x , 


) 


1.48391 x 10- 10 



Table 3: Approximation for Riccati matrix differential equation ( 16.13b in the interval [0,0.1] with step size h = 0.01. 



condition d6.9l l and taking h < 3/L, g(T) yields a contractive matrix function, which guarantees that equation (12.13| > 
has unique solutions A^. Hence, the matrix-cubic spline is completely determined. 

As an additional example for our proposed method, we consider the Riccati matrix differential equation ( 16. Il l with 

CW = { (l-x)x(2+x + 2x 2 ) l + (3-2x)x 2 +e*(x-x*) )> Y ^={ o)' (<U3) 

/ e x \ 

In this case, the problem has an exact solution given by Y(x) = I 2 ) > which will permit us to obtain the total 



x x 

error for all our numerical estimates. A short computation using expressions ( 16.4-b — (16.9b yields the following constants 



k Q = 6.13866 4o = 3 

w = V2 5 = 0.115758 

M= 12.0883 = 0.173205 

b = 2.23609 c= 1.17928 

d=\.Q\ L = 55.2443 



(6.14) 



which are necessary for the spline approximation in the interval [0,0.1], where 8 — 0.1 is taken for convenience. 
Therefore, we need to take h < 3/L = 0.0543042 and thus h = 0.01. The results are generated with Mathematica using 
FindRoot function to solve the emerging algebraic equations, and are summarized in Table 3, where the numerical 
estimates have been rounded to the fourth relevant digit. In each interval, we evaluated the difference between the 
estimates of our numerical approach and the exact solution, and then take the Frobenius norm of this difference. The 
maximum of these errors are indicated in the third column for each subinterval. 
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1.6e-10 




0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 



approx. interval 

Figure 3: Representing the absolute error margins for the Riccati matrix differential equation ( 16.131 ) in the interval 
[0,0.1] with step size h = 0.01. 

7 Conclusions 

This article develops a new method for the numerical integration of first-order matrix differential equations of the 
non-linear type Y'(x) = f(x,Y(x)),x £ [a,b] using matrix-cubic splines, and thereby generalizing the approach for the 
linear case in previous work (TJ. An important advantage of the proposed method is that the approximated solution is 
continuous in the interval under consideration, is easy to evaluate, and has an error of the order <9(/i 4 ). 
Our method is well-suited for implementation on numerical and/or symbolical computer systems (Mathematica, Mat- 
lab, etc.) as we have shown in Section 3 giving the explicit algorithm. For a full demonstration of our approach and its 
advantages, we conclude with two numerical examples for the Sylvester and Riccati matrix differential equations. 
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